Hello! Welcome to the preliminary exploratory analysis of using Rstudio for LiDAR data.
For those who don’t know. LiDAR data is known as light detection and ranging, it is a type of active remote sensing. Similar to other remotely sensed images, it can be collected via aerial, drone or satellite platforms.
In a nutshell, traditional maps you see are often passive remotely sensed data since it is easily accessible (i.e free, large-scale, variable locations and timescales, etc). However, LiDAR is becoming increasingly adopted as it collects information about earths surface across 3 axis (X, Y, Z), making it 3D.
FYI, passive remote sensing relies on sunlight and clear skies, active remote sensing relies on specific wavelengths (like microwaves) that can penetrate through clouds/ various weather conditions and isn’t impacted by night
Please refer to the image taken from SERC (2014).
Benefits of LiDAR that is increasing it’s adoption is that (and not limited to);
Issue regarding LiDAR is that it is not easily accessible due to it’s unique advantages. Thus, open source data is minimal and processing techniques are currently limited (i.e often on paid software). Which is why my department wanted to explore it’s potential processing via Rstudio.
My thesis and lab mates wanted to assess if the LiDAR data we collect from our in-lab made machines or open source files can be efficiently processed using free software to enhance our research. My thesis will be utilizing a ton of LiDAR data specifically in coastal ecosystems, which is why we wanted to assess the visualization and iteration capabilities of the LiDAR packages in Rstudio.
Packages for LiDAR data analysis are new developments and still with many bugs, missing functions, etc. The main package lidR was created in 2018 by Jean-Romain Roussel from Laval University. Till today he continuously works towards filling in the aforementioned gaps. The second specific package rLIDAR was created in 2021 by Silva, Crookston, Hudak, Vierling, Klauberg and Cardil from Moscow University.
A few tools and functions will be taken in this project to view the capabilities and limitations of using Rstudio for LiDAR data analysis. While trying to merge between the two.
To promote further push towards public data, data for this project is obtained from the minimal opensource LiDAR resources. Alberta is one of the few governments globally sharing their resources. Thus, data here is from Beaver Hills, Alberta.
to access their resources, kindly refer to this link: https://www.opendataareas.ca/#data
#Install Packages
#install.packages("lidR")
#install.packages("rLiDAR")
#install.packages("rgdal")
#install.packages("raster")
#install.packages("tmap")
#install.packages("tmaptools")
#install.packages("RStoolbox")
#install.packages("sf")
#install.packages("sp")
#install.packages("ggplot2")
#install.packages("rayshader")
#install.packages("mapview")
#install.packages("RMCC")
#install.packages("RCSF")
#Load packages
library("lidR")
library("rLiDAR")
library("rgdal")
library("raster")
library("tmap")
library("tmaptools")
library("RStoolbox")
library("sf")
library("sp")
library("ggplot2")
library("rayshader")
library("mapview")
library("RMCC")
library("RCSF")
las_beaver <- readLAS(paste0(getwd(),"/point_cloud_las/pc_083H07SW38.las"))
#view data
summary(las_beaver)
# check coordinate system
epsg(las_beaver) # checks coordinate system via points, if undefined will return zero
wkt(las_beaver) # checks coordinate system based on string rather than points, if undefined will be empty
# set coordinate system
st_crs(las_beaver) <-2955
las_beaver #check if it coordinate system is updated
# Determine if data is ready to be worked with, or if further wrangling required
las_check(las_beaver)
# some functions require the variable format to be LAScatalog, not Las file. So read in data in the two formats, for further processing.
# Read in file as LAScatalog (since there isn't a tool to convert it)
catalog_beaver <- readLAScatalog(paste0(getwd(),"/point_cloud_las/pc_083H07SW38.las"))
st_crs(catalog_beaver) <- 2956 #set coordinate system
catalog_beaver #check
plot(las_beaver) #las class (default is 3D)
plot(catalog_beaver, mapview = TRUE, map.type = "Esri.WorldImagery") #LAScatalog class(default 2D visualization)
The following is the las_beaver file (3D, interactive) please note the shape of the pixels
The following is the catalog_beaver file (2D and is overlayed on a map) please note the variation in file types
# Set Processing options
opt_chunk_size(catalog_beaver) <- 352 #352 x 352 meter tiles
plot(catalog_beaver, chunk_pattern = TRUE)
opt_chunk_buffer(catalog_beaver) <- 10 #overlap of 10 meter
plot(catalog_beaver, chunk_pattern = TRUE)
summary(catalog_beaver) #assess if processing extents set
# Create digital terrain model (DTM)
opt_output_files(catalog_beaver) <- paste0(getwd(),"/las_out/dtm_{XLEFT}_{YBOTTOM}") #store as raster
DTM <- grid_terrain(catalog_beaver, res = 4, knnidw(k = 10, p = 2), keep_lowest = FALSE, overwrite=TRUE) # cell-size set to 4x4 m, and k-nearest neighbor inverse distance weighting with 10 neighbors and a power of 2
# Visualize DTM
tm_shape(DTM)+
tm_raster(style= "cont", palette=get_brewer_pal("Greys", plot=FALSE))+
tm_layout(main.title = "Digital Terrain Model for Beaver Hill",legend.outside= TRUE)
writeRaster(DTM, paste0(getwd(),"/las_out/dtm.tif"))
#Histogram of Data Surface Elevation Distribution
DTM_distrubtion <- hist(DTM,
main = "Distribution of Surface Terrain Values (DTM) Pixels",
xlab = "Elevation (meters)", ylab = "Frequency",
col = "steelblue1")
# Create the components
slope <- terrain(DTM, opt='slope')
aspect <- terrain(DTM, opt='aspect')
hillshade <- hillShade(slope, aspect, angle=45, direction=180) #angle and direction based on the dataset being located in Beaver Hills, Alberta
# Visualize hillshade
tm_shape(hillshade)+
tm_raster(style= "cont", palette=get_brewer_pal("Greys", plot=FALSE))+
tm_layout(main.title = "Hillshade of Beaver Hills", main.title.position = "left", legend.title.size=0.1 , legend.text.size = 0.8, legend.position= c("center","bottom"), legend.outside = TRUE)
# Set processing exents for intensity analysis
opt_output_files(catalog_beaver) <- paste0(getwd(),"/lidar_points_out/intensity/int_{XLEFT}_{YBOTTOM}")
opt_filter(catalog_beaver) <- "-keep_first"
intensity <- grid_metrics(catalog_beaver, ~mean(Intensity), 5)
# Plot intensity
intensity[intensity<0]=0 # anything negative cannot be visualized and is equates to zero reflectance regardless
tm_shape(intensity)+
tm_raster(style= "quantile", n=5, palette=get_brewer_pal("Greys", n=5, plot=FALSE))+
tm_layout(main.title = "LiDAR Intensity in Beaver Hills", main.title.position = "left", legend.title.size=0.1 , legend.text.size = 0.8, legend.position= c("center","bottom"), legend.outside = TRUE)
Normalized digital surface model is used as it can separate between above ground and ground features. Often a prerequiste for canopy height model.
# Identify ground features
opt_output_files(catalog_beaver) <- paste0(getwd(),"/lidar_points_out/normalized/norm_{XLEFT}_{YBOTTOM}") #where files will be saved
normalized_beaver <- normalize_height(catalog_beaver, DTM) # spatial interpolation of features
# Identify above ground features (both built and natural) - Digital Surface Model
opt_output_files(catalog_beaver) <- paste0(getwd(),"/lidar_points_out/surface_model/dsm_{XLEFT}_{YBOTTOM}") #where files will be saved
DSM <- grid_canopy(catalog_beaver, res = 4, pitfree(c(0,2,5,10,15), c(0, 1))) #cell-size set to 4x4 m, and pitfree algorithim to smooth between features
writeRaster(DSM, paste0(getwd(),"/lidar_points_out/DSM.tif"), overwrite= TRUE)
# Calculate a Normalized Digital Surface Model
NDSM <- DSM - DTM
NDSM[NDSM<0]=0
NDSM
# Histogram of normalized digital surface model
NDSM_distrubtion <- hist(NDSM,
main = "Distribution of NDSM values",
xlab = "Elevation (meters)", ylab = "Frequency",
col = "coral")
# Visualize the Normalized Digital Surface Model
tm_shape(NDSM)+
tm_raster(style= "quantile", n=5, palette=get_brewer_pal("YlGn", n=5, plot=FALSE))+ # number of classes is 5 and using a quantile
tm_layout(main.title= "NDSM of Beaver Hill", main.title.position = "left", legend.title.size=0.1 , legend.text.size = 0.8, legend.position= c("center","bottom"), legend.outside = TRUE)
Looks specifically at treetops. Thus, not only the above ground features, but the highest point in the scene.
# Khosravipour et al. pitfree algorithm
thr <- c(0,2,5,10,15)
edg <- c(0, 1.5)
chm <- rasterize_canopy(las_beaver, 1, pitfree(thr, edg))
plot(chm)
# remove spurious effects of tree branches curtsey of Silva (2018)
data(chm) #changes format from spat raster to raster
chm_smooth <- CHMsmoothing(chm, "mean", 3, 0.67)
plot(chm_smooth)
This crashed the system but with the right server, when it ran, it gives you the individual classes/preliminary classification, which is highly beneficial. Should be carried out if possible.
# remove spurious effects of tree branches curtsey of Silva (2018)
data(chm) #changes format from spat raster to raster
chm_smooth <- CHMsmoothing(chm, "mean", 3, 0.67)
plot(chm_smooth)
Voxelization of pixels.
#3D voxel points
voxelize_p_beaver <- voxelize_points(las_beaver, 2)
plot(voxelize_p_beaver)
Notice that pixel shape has altered based on the feature, and the details provided (such as the lake bed). You can compare this to the image in the simple visualization tab.
Voxelization of objects rather than pixels.
# 3D voxel of objects
las1_dtm <- grid_terrain(las_beaver, res = 2, knnidw(k = 10, p = 2), keep_lowest = FALSE)
las1_n <- normalize_height(las_beaver, las1_dtm)
voxelize_beaver <- voxel_metrics(las1_n, ~sd(Z), res = 5)
voxelize_beaver
plot(voxelize_beaver, color = "V1", pal = heat.colors(60), mean) #visualization based on intensity in pixels
plot(voxelize_beaver, color = "Z", pal = heat.colors(60), mean) #visualization based on intensity in elevation
Horizontal transects shows variation along the X and Y axis.
#Create a horizontal transect
point1 <- c(371123, 5911101) #select first point
point2 <- c(371760, 5911101) # select second point
transect_h <- clip_transect(las_beaver, point1, point2, width = 5, xz = TRUE) #make the transect
# Visualize horizontal transect
ggplot(transect_h@data, aes(X,Z, color = Z)) +
geom_point(size = 0.5) +
coord_equal() +
theme_minimal() +
theme (legend.title=element_blank(),
axis.title.y=element_blank(),
axis.title.x=element_blank(),
axis.text.y=element_blank()
) +
scale_color_gradientn(colours = height.colors(50))
Vertical transects shows variation along, X, Y and Z axis.
# Create a vertical transect
point3 <- c(371270, 5911100) #select third point
point4<- c(371770, 5911477) #select fourth point
point5<- c(371789, 5912200) #check out variation with this point if interested
transect_v <- clip_transect(las_beaver, point3, point5, width = 2, xz = TRUE)
# Visualize vertical transect
ggplot(transect_v@data, aes(X,Z, color = Z)) +
geom_point(size = 0.2) +
coord_equal() +
theme_minimal() +
theme (legend.title=element_blank(),
axis.title.y=element_blank(),
axis.title.x=element_blank(),
axis.text.y=element_blank()
) +
scale_color_gradientn(colours = height.colors(100))
Another example;
#Set points for transect
pointA <- c(371123, 5911100)
pointB<- c(371789, 5912200)
#Filter ground points
ground_points <- filter_ground(las_beaver)
ground_transect <- clip_transect(ground_points, pointA, pointB, width = 4, xz = TRUE)
#Visualize ground points on transect
ggplot(ground_transect@data, aes(X,Z, color = Z)) +
geom_point(size = 0.5) +
coord_equal() +
theme_minimal() +
theme (legend.title=element_blank(),
axis.title.y=element_blank(),
axis.title.x=element_blank(),
axis.text.y=element_blank()
) +
scale_color_gradientn(colours = height.colors(50))
# All features on transect
ggplot(transect@data, aes(X,Z, color = Z)) +
geom_point(size = 0.5) +
coord_equal() +
theme_minimal() +
theme (legend.title=element_blank(),
axis.title.y=element_blank(),
axis.title.x=element_blank(),
axis.text.y=element_blank()
) +
scale_color_gradientn(colours = height.colors(50))
# Set function created by Jean-Romain Roussel (Laval University, Quebec, CA) as it is the current only way to combine transect location and plotting
plot_crossection <- function(las_class,
p1 = c(min(las_class@data$X), mean(las_class@data$Y)),
p2 = c(max(las_class@data$X), mean(las_class@data$Y)),
width = 2, colour_by = NULL)
{
colour_by <- enquo(colour_by)
data_clip <- clip_transect(las_class, p1, p2, width)
p <- ggplot(data_clip@data, aes(X,Z)) + geom_point(size = .7) + coord_equal() + theme_minimal()
if (!is.null(colour_by))
p <- p + aes(color = !!colour_by) + labs(color = "")
return(p)
}
# Using Progressive Morphological Filter (PMF) created by Zhang et al. (2003)
las_class_A <- classify_ground(las_beaver, algorithm = pmf(ws = 4, th = 2))
plot_crossection(las_class_A, point3 , point4, colour_by = factor(Classification))
# Using Multiscale Curvature Classification (MCC) created by Evans & Hudak (2017)
las_class_B <- classify_ground(las_beaver, algorithm = mcc(s = 4, t = 0.3))
plot_crossection(las_class_B, point3 , point4, colour_by = factor(Classification))
# Using Cloth Simulation Filter (CSF) created by Zhang et al. (2016)
las_class_C <- classify_ground(las_beaver, algorithm = csf(TRUE, 0.5, 1, time_step = 0.65))
plot_crossection(las_class_C, point3 , point4, colour_by = factor(Classification))
Please read the research paper. To access files and R processing, github repository can be access here: https://github.com/lanajbb/ES8913_lidR_project
Gomez, C. (2020, December 30). Lidar data with the lidR package - 1. https://rstudio-pubs-static.s3.amazonaws.com/708624_40f54f88124b4b6eb214006a3a53a691.html
Luo, W., & Zhang, H. (2015). Visual analysis of large-scale LiDAR point clouds. Proceedings - 2015 IEEE International Conference on Big Data, IEEE Big Data 2015, 2487–2492. https://doi.org/10.1109/BIGDATA.2015.7364044
Roussel, J. R., Auty, D., Coops, N. C., Tompalski, P., Goodbody, T. R. H., Meador, A. S., Bourdon, J. F., de Boissieu, F., & Achim, A. (2020). lidR: An R package for analysis of Airborne Laser Scanning (ALS) data. Remote Sensing of Environment, 251. https://doi.org/10.1016/J.RSE.2020.112061
SERC. (2014, July 15). Figure 2. Satellite Imagery versus LiDAR. https://serc.carleton.edu/details/images/47158.html
Silva, C. A., Hudak, A. T., Vierling, L. A., Valbuena, R., Cardil, A., Mohan, M., Almeida, D. R. A., Broadbent, E. N., Almeyda Zambrano, A. M., Wilkinson, B., Sharma, A., Drake, J. B., Medley, P. B., Vogel, J. G., Prata, G. A., Atkins, J. W., Hamamura, C., Johnson, D. J., & Klauberg, C. (2022). treetop : A Shiny‐based application and R package for extracting forest information from LiDAR data for ecologists and conservationists . Methods in Ecology and Evolution. https://doi.org/10.1111/2041-210X.13830
Silva, C. A., Klauberg, C., & Mohan, M. M. (2018). LiDAR Analysis in R and rLiDAR for Forestry Applications rGEDI: An R Package for NASA’s Global Ecosystem Dynamics Investigation (GEDI) Data Visualization and Processing. View project GEOAGRO-Geoprocessamento e Sensoriamento Remoto Aplicados a Sistemas Agroecológicos e Mapeamento da Cobertura Vegetal e Uso da Terra (MCVUT) View project. https://www.researchgate.net/publication/324437694
Wasser, L., & Jones, M. (2021, May 13). Data Activity: Visualize Elevation Change using LiDAR in R to Better Understand | NSF NEON | Open Data to Understand our Ecosystems. https://www.neonscience.org/resources/learning-hub/tutorials/da-viz-neon-lidar-co13flood-r
Zhou, T., & Popescu, S. (2019). remote sensing waveformlidar: An R Package for Waveform LiDAR Processing and Analysis. MDPI. https://doi.org/10.3390/rs11212552